#!/usr/bin/Rscript

library(fields)
h0 = 0 #bottom elevation above datum
Sf = 10 #friction slope
Sy = 1 #specific yield
K =  .1 #wetland conductivity (m/s)
alpha=.5; beta=4 #kadlec coefficients
L = 100 #length of pond (m)
w = 12 #width of pond (m)
nodes = 1000 #number of grid points
dx = L/(nodes-1) #spatial discretization (m)
dt = 1 #temporal discretization (sec)
num.t = 3600 #number of time steps

#initial and boundary conditions
h  = array(1,dim=c(nodes,num.t)) #initial head above datum (m)
h[1,] = array(1,dim=num.t) #left boundary in time
h[nodes,] = array(.5,dim=num.t) #right boundary in time

#simulation
my.colors=array(tim.colors(),dim=num.t)
plot.new()
plot.window(c(1,nodes),c(0,2))
axis(1,seq(from=1,to=nodes,by=10))
axis(2,c(0,1,2))
lines(x=(1:nodes),y=h[,1],col=my.colors[1])
for (t in 2:num.t) {
	for (i in 2:(nodes-1)) {	
		dummy = (h[i+1,t-1]-h[i,t-1])*(h[i+1,t-1]-h0)^beta -  (h[i,t-1]-h[i-1,t-1])*(h[i-1,t-1]-h0)^beta
		h[i,t] = dummy*(K*Sf^(alpha-1)*dt)/(2*dx*Sy)+h[i,t-1]
	}
	lines(x=(1:nodes),y=h[,t],col=my.colors[t])
	print(t)
}

	




 
